Ground state properties of the one dimensional Coulomb gas. 
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We study the ground state properties of a quasi one dimensional electron gas, interacting via 
an effective potential with a harmonic transversal confinement and long range Coulomb tail. The 
exact correlation energy has been calculated for a wide range of electron densities by using the 
lattice regularized diffusion Monte Carlo method, which is a recent development of the standard 
projection Monte Carlo technique. In this case it is particularly useful as it allows to sample 
the exact ground state of the system, even in the low density regime when the exchange between 
electrons is extremely small. For different values of the width parameter 6 (0.1 aJ5 < b < 4 aj$), we 
give a simple parametrization of the correlation energy, which provides an accurate local density 
energy functional for quasi one dimensional systems. Moreover we show that static correlations 
are in qualitative agreement with those obtained for the Luttinger liquid model with long range 
interactions. 



I. INTRODUCTION 

The recent experimental realizations of ideally clean 
quasi one dimensional (Q1D) systems, like ultra cold 
Fermi gases in elongated harmonic traps^ and high mo- 
bility quantum wires in the so-called cleaved edge over- 
growth samples 2 -, have stimulated a new intense theo- 
retical effort to explain the physical outcome of these 
systems. The low dimensionality brings about pecu- 
liar phenomena such as the fractionalized (0.7 struc- 
ture) conductance 3 -^, enhances the effect of localization 
(Wigner crystallization)^, breaks the validity of the 
Fermi-liquid paradigm, which must be abandoned in fa- 
vor of the Tomonaga-Luttinger liquid (TLL) concept 6 . 
One of the most striking consequences of the Luttinger 
theory is the spin-charge separation, which has been seen 
in a series of remarkable experiments carried out by Aus- 
laender et ai, who were able to resolve the dispersion en- 
ergy of elementary spin and charge excitations^ 7 -. They 
used the tunneling current between two parallel wires 
to probe the properties within one of the two wires. 
However, their measurements exhibit some features, like 
fringes in the tunneling pattern and non unitary conduc- 
tance, that are not completely understood. 

The finite-size effects, the disorder and the inhomo- 
geneity of the device can play a crucial role to quanti- 
tatively explain the experiments 8 -^^. In this paper we 
will rather focus on the simpler homogeneous Q1D elec- 
tron gas with harmonic transversal confinement and ef- 
fective interactions with a long range Coulomb tail (1/r). 
The details of the confinement only affect the behav- 
ior at short range of the effective potential, and many 
model o 11 ' 12 ! 13 ! 14 ! 15 have been proposed which give an 
equivalent description of the homogeneous Q1D Coulomb 
gas. Despite its simplicity and a huge amount of theoret- 
ical wor k 16 i 17 ' 18 i 19 i 20 done to understand its properties, 
an accurate parametrization of its correlation energy is 
still missing. Indeed, the calculation of the ground state 



energy of ID wires with realistic Coulomb interactions is 
still an open problem, since the TLL is an effective low 
energy theory, and the RPA perturbative expansion is 
correct only in the high density limit. Recently a map- 
ping of the problem with realistic Coulomb interaction 
onto exactly solvable models has been proposed 2 ^, but 
within this scheme several approximations are required 
for different density regimes. Previously an STLS-like 23 
method was used by Calmels and Gol d 22 i 24 to compute 
the correlation energy, but it turns out to be not accurate 
enough to yield the correct ground state in the low den- 
sity regime. Indeed, it predicts a Bloch instability ruled 
out by the Lieb-Mattis theorem 2 ^. 

Projection quantum Monte Carlo (QMC) techniques 
are exact in one dimension, since in this case the ground 
state (GS) nodes are known exactly^ 6 -, and the so called 
fixed node approximation, which cures the well known 
sign problem, does not affect the results. However, pre- 
vious diffusion Monte Carlo (DMC) simulations^ 7 - suf- 
fered from a lack of ergodicity at low electron densities, 
when the exchange between electrons become exponen- 
tially small. Other QMC attempts^ to study the Q1D 
electron gas used the "world-line" algorithm on the lat- 
tice after a naive discretization of the Laplacian. Here we 
apply the novel lattice regularized diffusion Monte Carlo 
(LRDMC)^ 9 ., which is more efficient than previous pro- 
jection QMC methods at low density, where it substan- 
tially alleviates the lack of exchanges between electrons. 

The aim of this work is to provide a simple and ef- 
ficient parametrization of the ground state correlation 
energy, exactly computed for the unpolarized system. 
Other ground state properties are also studied, like the 
spin and charge structure factors, which reveal a strong 
similarity with those computed using a TLL with long 
range interactions^ 6 -. Therefore our parametrization can 
be an extremely useful input of density functional theory 
(DFT) calculations of Q1D systems with local density 
approximation (LDA), which can include the homoge- 
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neous Q1D electron gas with realistic long range interac- 
tions and TLL features as the reference system. Previ- 
ous successful attempts have been limited so far to model 
systems such as the Luttinger liquid 3 ^ or the one dimen- 
sional Hubbard model^L, that have been successfully used 
as the reference systems for DFT simulations of Q1D ul- 
tra cold inhomogeneous atomic gases. We believe that 
an essentially exact calculation of the correlation energy, 
presented in this paper, should open the way for a wide 
range of relevant realistic applications in the field of Q1D 
systems. 

The paper is organized as follows. In Chap. [IT] we 
present the model for the Q1D homogeneous electron 
gas, in Chap. [TTTI we describe the variational ansatz used 
in our simulations, in Chap. IIVI we briefly review the 
LRDMC method and we compare it with the standard 
DMC algorithm. The results are reported in Chap. [Vj 
where we present the parametrization of the correlation 
energy with different values of the transversal confine- 
ment, and in Chap. IVI1 where we study the charge and 
spin structure factors. Finally, the conclusions are drawn 
in Chap. IVIIl In the end, two appendices explain how to 
compute the RPA correlation energy at high densities, 
and how to estimate the plasmon excitations from the 
knowledge of the static structure factor. 



II. MODEL 

In this paper we study a realistic model for a quan- 
tum wire with the lateral confinement provided by a har- 

2 

monic transversal potential V(r±) = j^, where b tunes 
the strength of the confinement and measures the wire 
width. Here and henceforth we use the effective Bohr 
radius ajj = " e a as length unit and the effective Ryd- 

2 

berg Ryd* = as energy unit, where e is the dielectric 
constant of the semiconducting medium and m* is the ef- 
fective mass of the electrons in the semiconductor. The 
electrons in the wire interact via a long range Coulomb 
potential. If the confinement is sufficiently strong, the 
ground state (GS) of this system can be approximated 
with good accuracy by a wave function with longitudinal 
and transversal components factorized. In particular we 
neglect any contribution from higher subbands of the lat- 
eral direction and we take the GS of the two dimensional 
harmonic oscillator as the transversal part of the total 
wave function. This approximation is valid whenever 
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where r s is the Wigner-Seitz radius (2r s = 1/p is the 
mean interparticle distance), i.e. for sufficiently low elec- 
tronic density. 

Tracing out the transverse motion from the full 
Schrodinger equation by integration over the lateral co- 
ordinates of the particles yields^, for N particle on a 
segment of length L, the one dimensional (ID) Hamilto- 
nian 



N 1 

H = E V * + 2L £ Vb ^ - ^ 



(2) 



with p(k) — ^2 ■ exjp(ikxj) the Fourier transform of the 
one-body density operator and 



V b (k) = 2E 1 {b 2 k 2 )eMb 2 k 2 ). 



(3) 



Above E\ is the exponential integral function and in Eq. 
((2|) a suitable rigid positive charge background exactly 
cancels the k = term. The real space ID interparticle 
potential, 



(4) 



has a long-range Coulomb tail but is finite at the origin. 

Since in this work we focus on the GS properties of 
the wire, we work in the sector of vanishing total spin 
component along the z quantization axis, namely = 
= N/2. Obviously this choice is not a restriction 
because the ground state certainly belongs to this sector, 
regardless of its total spin. On the other hand the value 
of the total spin is known from the Lieb-Mattis theorem 
and Ref. [2?]: the GS of the one dimensional Coulomb gas 
is always a singlet for all densities. In order to perform 
a finite-size scaling analysis of the energy and correlation 
functions, we carried out QMC simulations with different 
number of particles, going from N — 10 to N = 162 and 
with an odd number of particles per each spin so that the 
degeneracy effects are avoided. With the aim to further 
reduce the finite-size bias, we considered the Hamiltonian 
in Eq. [2] with periodic boundary conditions (PBC) and 
with an infinite number of replicas of the simulation box 
(supercell). Thus, an electron in a supercell interacts 
with the other electrons in the supercell, their images, its 
own images, and the background. It is then convenient 
to define an effective interparticle potential, by summing 
the bare interaction of a particle and its background with 
a second particle over all its images, to obtain a periodic 
function: 



V(x) = 



V b (x + nL) - 



L/2 



L/2 



dy Vb(x + nL-y) 



(5) 



In the above expression L is the length of the simulation 
box, n takes relative integer values and G n — 2i:n/L is a 
reciprocal vector of the ID Bravais lattice with primitive 
unit cell of length L. Since Vb is a long-range potential, 
we resort to an Ewald-likc method 32 to compute the sum 
in Eq. [5l The short-range part of the potential V and 
its long-range tail are treated in a different fashion, the 
former being summed in the direct space, the latter in 
the reciprocal one. The Ewald's procedure yields: 



3 



V(x) 
V sr (x) 



V SI (x) + V lr (x), 



—j^- ^ exp (x — nL) 2 /4b 2 erfc 

n— — oo 



nil 



2b 



+ 00 

E 



2 J\x-nL\ 
7 en 1 



\x — nL\ 



2b 



(6) 
(7) 

(8) 



Tl>0 



In practice, we have worked with the Hamiltonian 

N N 

^ = -E V ' + E^) + Y^MAD, (9) 
i=l i<j 

where Vmab = V(0) — V&(0) is the Madelung energy, i.e., 
the interaction of a particle with its own images. We have 
used a tabulation for the potential V(x). In particular, 
the G sum in Vi T has been truncated at G = 12/6 and 
a sufficient number of images in V sr has been included, 
so that the overall error in the tabulation is less than 
1Q- 6 Ryd*. 

III. WAVE FUNCTION 

The wave function we used in our QMC simulations 
is of the Slater-Jastrow type: 

*t - exp f-5^«(a:y)j # 1 (10) 

where D a is a determinant of N a plane waves with wave 
vectors occupied up to the Fermi momentum kp = 
The Jastrow factor contributes significantly to improve 
the quality of the variational state, since it correlates the 
particles and tunes the amplitude of ^t- Here we use 
a two body Jastrow factor, which takes into account the 
electron-electron correlation, without spoiling the trans- 
lational invariance of the system. A recent work 3 ^ on the 
ID t — t' Hubbard model has shown that the long-range 
behavior of the two body Jastrow effectively accounts 
for the proper description of the metallic and insulating 
phases of that lattice model. Therefore a good functional 
form of u(x) in Eq. [TO] is crucial to obtain the correct 
physics for a strongly correlated ID system. 

In order to avoid spin contamination 3 —, the function 
u{x) does not depend on the spin of the particles. In 
particular, we choose the RPA form of u(x) as our first 
variational ansatz. Following Ref. [3f| the Fourier com- 
ponents of u are: 

2pu RPA (k) =-S (k)- 1 + VSo(k)- 2 + 2pV b (k)/ki, 

(11) 

with S (k) = (k/2k F )9(2k F - k) + 9(k - 2k F ) the struc- 
ture factor of a non interacting one dimensional electron 



gas (1DEG). The two body Jastrow is repulsive at the 
origin, to reduce the probability of two electrons to ap- 
proach each other, and lower in this way the potential 
energy. Since the effective Coulomb coupling increases 
as the density decreases, the lower is the density, the 
higher is the Jastrow repulsion. 

Rescaling urpa is the most straightforward improve- 
ment beyond the RPA ansatz. It improves the variational 
wave function by means of the simple parametrization 

u(x) = ju RPA (x), (12) 

where 7 is a linear parameter that we optimized using 
the variance minimizatio n 3 ^ 37 . The optimal values for 7 
are reported in Tab. Q] for different densities and for N = 
22: the RPA seems good for the highest densities, where 
7 ~ 1, but it becomes worse for lower densities, where the 
rescaling is effective. Indeed, whenever the correlation is 
stronger, urpa overestimates the interelectron repulsion. 
Moreover we have found that the optimal value of 7 does 
not depend on the number of particles in the supercell, 
the function urpa having the proper dependence on N. 



TABLE I: Scaling parameter 7 of the urpa function op- 
timized using the method of variance minimization for the 
1DEG with b = 0.1 and N = 22 
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With the aim to check whether the scaled RPA is ac- 
curate enough to yield the correct Jastrow correlation 
for the lowest densities, when the effective coupling is 
higher, we take into account the most general expression 
for the two body Jastrow. We expand u in a linear sum 
of Chebyshev polynomials^, which are a complete basis 
set in the orthogonality interval (—1,1): 

u ( x ) = E a ^ T 2m I ~ I 2 j , (13) 
m>0 \ 2 / 

where the range (0, j) is mapped into the interval (—1, 0) 
and only the even polynomials are used. In this way the 
condition u'(L/2) = is fulfilled, and ensures the conti- 
nuity of the first derivative of u at the edge of the super- 
cell. Moreover the sum in Eq.[T3] starts from m = 1, since 
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To is the identity and the wave function is determined 
apart from a constant. Using the variance minimization, 
we optimize the parameters a m up to the convergence 
of the expansion in Eq. 1131 which has been reached for 
m = 10. The variational energies relative to the vari- 
ous functional forms of u for r s — 10 and N — 22 are 
summarized in Tab. [TT1 whereas the functions are plotted 
in Fig. [T] The comparison shows that the scaled RPA 
Jastrow factor leads to a very good variational state, as 
its energy is very close to the exact GS value (more than 
99.5% of correlation energy is recovered) and it is almost 
coincident to the most flexible variational form obtained 
with the mentioned Chebyshev expansion. The simple 
RPA form is a good approximation, since it provides a 
large fraction of the correlation energy (98.3%), but the 
rescaling of the RPA Jastrow yields a further substan- 
tial improvement. Thus it is clear that the most conve- 
nient parametrization of the Jastrow term is the scaled 
urpai because one is able to reach an almost exact wave 
function already at the variational level, by optimizing 
just one variational parameter. For the above reason, in 
the forthcoming sections we will use the optimized and 
rescaled RPA Jastrow as trial wave function ^>t for all 
densities (r s ) and sizes {N) taken into account in this 
QMC study of the 1DEG. 
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FIG. 1: Optimized u functions for r s — 10 and N = 22: 
urpa{x) (dotted line), ^urpa{x) (dashed line), and the 
Chebyshev expansion for u(x) (solid line). r„a,Q is the unit 
length. 



IV. LATTICE REGULARIZED DIFFUSION 
MONTE CARLO METHOD 

The low dimensionality and the strong correlation 
among the electrons not only have dramatic consequences 
on the physical properties of the quantum wire, which 
will be studied in Sees. [V] and IVI1 but also affect the 
efficiency of the QMC simulations of the system. In- 
deed, as we have seen in the previous section, the effective 



TABLE II: Total energy E to t, correlation energy E C orr and 
percentage of correlation energy %E corr for b — 0.1, r s = 10, 
and N — 22. The fraction of the correlation energy recov- 
ered is computed from LRDMC calculations which provide 
the exact GS energy for a 1DEG (see Sec. HVfl . 





Etot 


E 

J-^corr 




RPA 


-0.47207(2) 


-0.20519(56) 


0.9830(15) 


Scaled RPA 


-0.474825(9) 


-0.20794(55) 


0.9962(15) 


Chebyschev 


-0.474900(9) 


-0.20802(55) 


0.9965(15) 



Coulomb interaction leads to an optimal wave function 
with a strong repulsive Jastrow factor, which freezes the 
relative positions of the particles and enhances the 2kp 
component of the charge-charge correlation, a signature 
that the system is close to the Wigner phase. The two 
body Jastrow, necessary to provide a good variational 
description of the system, introduces pseudo nodes, i.e. 
surfaces in the configuration space where the wave func- 
tion almost vanishes, due to the exponentially increase 
of J(x), which acts like a Gutzwiller projector, by avoid- 
ing "double" occupancies on a given electronic position 
x. These pseudo nodes are similar to the usual nodes of 
the fermionic wave function, but if the latter arise from 
the antisymmetrization of the many body state, the for- 
mer are a consequence of the strong repulsion, which pre- 
vents two electrons to come closer and eventually overlap. 
The effect of the pseudo nodes on the QMC simulation is 
harmful, since in the ID case they can lead to a slow con- 
vergence of the Markov chain to the equilibrium distribu- 
tion. In particular, it can be extremely difficult to con- 
nect two configurations with a spin exchange. However, 
the charge degrees of freedom are still well reproduced in 
spite of this lack of ergodicity, as one might infer that the 
pseudo nodal pockets are equivalent for the charge prop- 
erties. Instead the expectation values of spin dependent 
operators are spoiled, if the Markov chain is not able to 
guarantee a sufficient number of spin exchanges during 
the simulation in a reasonable time. 

The variational Monte Carlo (VMC) algorithm can 
easily overcome the problem, since the proposed move 
can be forced to flip the spin of an electronic configura- 
tion, either by explicitly introducing a spin exchange or 
by allowing the amplitude of the move to be greater than 
the mean interparticle distance. Instead in the diffusion 
Monte Carlo (DMC) approach, the random walk has to 
follow the diffusion process driven by the imaginary time 
dependent Schrodinger equation. If the importance sam- 
pling is introduced, the resulting Green function, approx- 
imated by means of the Trotter expansion up to the first 
order in the time step r, includes the drift-diffusion dy- 
namics: 

R' = R + L>TVln|* T (R)| 2 + X V2Dt, (14) 

where V In ^^(R)! 2 is the quantum force, D — 1 is the 
diffusion coefficient, and x is a Gaussian distributed ran- 
dom variable. The configurations generated step by step 
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are distributed accordingly to \^>t\ 2 at the beginning of 
the simulation, but after a transient they will reach the 
equilibrium and sample the mixed distribution ^fn^t, 
where VP fn is the lowest variational state with the same 
nodes as the trial wave function (this is the so called 
fixed node (FN) constraint). In order to get rid of the 
time step bias in the final result, one needs to extrap- 
olate the FN energies obtained at different time steps 
for t going to zero. For smaller r, spin exchanges be- 
come rarer, because the mean square displacement van- 
ishes linearly with r (Eq. [T4|) . To overcome the lack of 
spin exchanges in the DMC, we have applied a different 
projection QMC method, the lattice regularized Green 
function Monte Carlo (LRDMC), successfully introduced 
in Ref. to cure the localization error in the presence 
of non local potentials. In this section we review the 
method and compare its efficiency of sampling spin flips 
with respect to the standard DMC framework. 

The main idea behind the LRDMC method is to 
deal with a regularized Hamiltonian in such a way 
that the standard Green function Monte Carlo (GFMC) 
algorithm 3 — i 40 ' 41 on a lattice can be applied also to con- 
tinuous systems. The regularization of the Hamiltonian 
in Eq. [5] involves both the kinetic and the potential parts. 
The Laplacian is discretized by means of the finite differ- 
ences 

A = i] \ P A a + (1 - p)A a '] +0(a 2 ), (15) 

with A a an Hermitian lattice operator given by 

A a ^( Xl ) = \ (V{ Xi + a) + *( Xi -a)- 2V{ Xi )) , (16) 

where a is the mesh size, p and r\ are constants (r\ — 
1 + 0(a 2 )), and Xi is the position of the i-th electron. 
Due to the homogeneity of the system, p is kept spatially 
independent, contrary to the general case^ where the 
dependence of p on the electronic positions is exploited 
to improve the efficiency of the diffusion process. Here 
p = 0.5, and the contributions to the total Laplacian 
coming from A a and A a are equally weighted. The two 
terms, with a' /a = allow the diffusion to explore all 
the continuous space, since the two meshes are incom- 
mensurate; in this way the lattice space bias due to the 
discretization of the continuous kinetic operator is greatly 
reduced and one can work with a reasonably large value 
of a without a significant lattice step error. 

Also the potential is regularized, so that our final 
Hamiltonian H a fulfills the following three conditions: 
i) H a — > H for a — > 0; ii) for the chosen trial wave func- 
tion ^t, for any a and any configuration x, the local 
energy ei(R, [^t]) = ^J/r" of the continuous Hamilto- 
nian H is equal to the one e£ (R, [^t] ) corresponding to 
the Hamiltonian H a ; iii) the discretized kinetic energy is 
equal to the continuous one calculated on the state ^t- 
The condition (iii) determines the constant 77, while the 
condition (ii) fixes the form of the regularized potential 



ya. 



V a (R) = V(R) 



(R) 



(17) 



Notice that the condition (ii) yields another important 
property for H a : if is an eigenstate of H, it is also 
an eigenstate of H a for any a. Thus, as the quality of 
^>t increases, the dependence of the LRDMC energy on 
a decreases. 

The lattice regularized Hamiltonian H a reads: 



^R'.R - 



—77 p/a 2 

-r) (1 -p)/a' 2 



if R = R + 5 a 
if R = R + S a > 
V a (R) ifR = R, 

(18) 

where 5 a (5 a ') is a N dimensional vector defined as the 
one particle displacement of length ±a (±a'). Thus 
there are 2N different 5 a (S a > ), and H a in Eq. fTS] con- 
tains AN off diagonal elements, which come from the 
discretization of the Laplacian (Eq. I15p . In particu- 
lar, by defining the importance sampling Green function, 
Gr,, r = <Jt(R)(A<W - #r',r)/*t(R), the config- 
uration R is connected by Gr'.r to a finite number of 
configurations R', although R and R' live in a contin- 
uous space. Therefore the Green function Gr^r is dis- 
crete, and can be sampled using a heat bath algorithm, 
like in the standard GFMC scheme on a lattice, although 
in this case R and R' are continuous variables. Another 
important difference with respect to the lattice case is the 
spectrum of H a , which is unbounded from above; thus, 
in order to guarantee the positivity of the Green func- 
tion, we need to perform the limit A — * 00, which can be 
handled within the continuous time formulation, already 
introduced in Ref. EH for the GFMC method. 

Although in the continuous limit a — » there is no 
sign problem because the sampling is restricted within a 
region -the nodal pocket- with definite sign, for non zero 
a the fermionic sign problem is still present and needs 
to be treated by means of an effective Hamiltonian 39 , 
which approximates the regularized H a (FN constraint). 
The regularized effective Hamiltonian H cS included in 
the LRDMC algorithm is defined as follows: 



rreff 
H R,R 



Tja 

H R,R> 







Tja 

W R,R 



V./(R) 



if R ^ R' and 
*t(R)#rr'/*t(R) < 
if R ^ R' and 
$ T (R')^ R ,/f T (R) > 
if R = R', 

(19) 

where V s/ (R) = Er ¥ r*t(R')^,r/*t(R) > 0, the 
so called sign-flip term, is the sum over all the terms 
that cause a negative sign problem in the Monte Carlo 
sampling. For this reason these terms are traced in the 
diagonal part of the effective Hamiltonian, that therefore 
no longer contains off diagonal terms with the "wrong" 
sign. The ground state ^fn of H cS has the same signs 
as the trial wave function Wt, and so the mixed distri- 
bution 'i'FN'&T sampled during the LRDMC simulation 
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will be non negative. In general, in the limit a — > 
the FN energy Efn — {^fn\H cS \^t) of the effective 
Hamiltonian H cS is an upper bound of the GS energy Eq 
of H. The FN approximation turns out to be exact only 
if the nodes of the trial wave function are the same as 
the GS nodes. However, the trial wave function in 
Eq. [TU] has the exact GS nodes. Indeed in the ID case 
the nodal structure of the GS is exactly defined by the 



coalescence planes Xi 



where Xi and two elec- 



trons with the same spin, and the position of the planes 
are completely determined by the antisymmetry of the 
particles 26 . Therefore both the LRDMC (a -> 0) and 
DMC (r — * 0) values for this ID system are exact within 
their statistical precision. 

We did an accurate comparison between the DMC 
and LRDMC approaches, by taking into account the ef- 
ficiency of the energy estimate, the dependence on the 
time step and on the lattice space, and the spin flip fre- 
quency during the simulations, defined as the number 
of exchanges between two particles with opposite spin 
per unit time (imaginary projection time) per particle. 
The "standard" DMC algorithm, taken as the reference 
for our comparison, is described in Ref. We applied 
the two QMC schemes to the quantum wire model with 
N = 22 and for r s ranging from 1 to 10. For a fair com- 
parison, we chose the DMC time step r = a 2 /2, so that 
both algorithms provide the same amplitude for the diffu- 
sion move. The efficiency of the DMC energy estimate is 
double that of the LRDMC, since in the latter approach 
we need to compute in advance all the possible off diago- 
nal moves, thus losing a fraction of the computing time. 
On the other hand, as reported in Table IIII1 the spin ex- 
change frequency is almost the same for the high density 
model, when the correlation is weak, but the LRDMC be- 
comes more and more effective in sampling the spin flips 
when the density is lower and the correlation becomes 
stronger. In particular, for the lowest density (r s = 10) 
the LRDMC algorithm yields an efficiency in the spin flip 
sampling which is two orders of magnitude higher than 
the DMC case. 

Another appealing behavior of the LRDMC approach 
is the lattice space dependence of the fixed node energy. 
As one can see in Fig. [2 the LRDMC energies have a 
quadratic dependence on a with a prefactor much smaller 
than the slope of the linear fit for the corresponding DMC 
energies. This means that in order to obtain an almost 
converged LRDMC result one does not need to go to 
small lattice spaces, with a gain both from the compu- 
tational point of view and from the efficiency in the spin 
flip sampling, which of course is reduced as the diffusion 
move goes to zero. 

To summarize, it is apparent that both the lattice step 
bias and the lack of ergodicity are greatly reduced by us- 
ing the LRDMC algorithm in the place of the standard 
DMC. Therefore, given the amplitude of the QMC move, 
the LRDMC is more effective than the DMC scheme. We 
believe that the reason is related to the Trotter approx- 
imation behind the DMC propagator, which spoils the 
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For both the two cases, the dependence appears to be linear, 
with a slope of 0.117(8) for the DMC algorithm and -0.018(2) 
for the LRDMC approach. 



exact dynamics of the diffusion process and apparently 
affects the ergodicity of the random walk. On the other 
hand, even within a finite lattice space, the LRDMC algo- 
rithm converges to the exact GS of the effective Hamilto- 
nian H in Eq. 1101 implying that simulations are always 
physically meaningful even for a > 0. This should allow 
in general a more controlled and smooth extrapolation to 
the continuous a — > limit. 

In the next section, we will analyze both the VMC and 
the LRDMC results for the quantum wire model. For ev- 
ery density and width, we have performed a lattice space 
extrapolation (a — > 0) and a finite size extrapolation to 
the thermodynamic limit (N — > oo). To extrapolate the 
energies to the continuous limit, we have fitted points 
computed in the range 0.05 < a < 1.2 with a quadratic 
function in a, as reported in Fig. [51 while for the extrap- 
olation to the thermodynamic limit we used the function 
E OQ +c 1 /N+c 2 /N 2 in our fits. We evaluated the LRDMC 
energies for N — 10, 22, 42, 62, 82, but in some cases, for 
low r s , we carried out LRDMC simulations with as many 
as 242 particles, in order to have always a reliable esti- 
mate of the finite size errors. 



V. CORRELATION ENERGY 

The correlation energy per particle E corr = Eq — Ehf 
is computed for various values of the width parameter b 
and densities r s , and parametrized as a function of r s for 
each b. The unpolarized HF energy of this quantum wire 
model is 



E HF (r s ,b) 



48r? 



(20) 
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TABLE III: Spin exchange frequency (Ryd*) for the LRDMC and DMC algorithm at different densities for the quantum wire 
model with N — 22 and b — 0.1. Notice that the frequency is reduced when the density decreases, while the efficiency of the 
LRDMC increases with respect to the DMC. All the simulations have been performed with a = 0.2r s and r = a 2 / 2. 





r s = 1 


r s = 2 


r s = 4 


r a = 6 


r s = 10 


LRDMC 


1.18 


5.72 10~' J 


1.96 lO - ^ 


3.19 10~ 4 


8.86 10~ b 


DMC 


1.14 


3.81 10^ 


4.24 10~ 3 


3.21 10" b 


9.09 10~ B 


rel eff 


1.03 


1.50 


4.62 


9.94 


97.47 



where the first term is the kinetic energy, and the second 
contribution is the exchange energy (E ex ) term with a = 
4/ir and the function F defined by 



F(R) = 



1 

2^ 



l/R 



dxf(x) [1 - Rx] . 



(21) 



The Hartree term vanishes since the system is neutral 
and homogeneous. 

We evaluated the GS energy E using the LRDMC 
method with the FN approximation, which projects the 
initial to the lowest energy state of the system with 
the same nodes of ^t- However, since the trial wave 
function in Eq. [TO] has the exact GS nodes, our LRDMC 
energies are exact within their statistical precision in the 
limit a — > 0, as already pointed out in Sec. IIV1 The 
correlation energy is then obtained by subtracting the 
HF energy of Eq. [20] from the GS energy. Thus, the 
correlation energy is computed exactly for a given value 
of r s and b. 

It is useful to study the correlation energy in the 
high and low density limits, in order to find a good 
parametrization which includes the correct asymptotic 
behavior. In the high density limit, i.e. r s — > 0, the 
correlation energy can be computed via a perturbative 
expansion of the interaction, using the RPA technique 
to find the coefficient of the lowest order term in r s (see 
Appendix [K§ . It turns out that the correlation energy is 
quadratically vanishing, as r s goes to zero 



r(r s -> 0) = 



.4 



corr 9 
T 

• a 1 



ir 4 b 2 



(22) 



where A corr = J °° dx xf(x) 2 — 4.9348. This result was 
obtained by Calmels and Gold, using the mean spherical 
approximation (MSA) 22 ' 23 , which is consistent with the 
RPA finding. On the other hand, in the low density 
regime (r s — ► oo) the exact behavior of E corr can be 
guessed by studying the ratio E corr / E ex . For instance, 
in Ref. [22| the correlation was computed using a three- 
sum-rule approach (3SRA) of the STLS theory for the 
same model, and this ratio turned out to be 



Ecorr ys * oo) 



0.84. 



(23) 



E ex {r s -* oo) 

Since the asymptotic expansion of E ex (r s ) in the limit of 



£tl( , ) ^Q + 0( _l_ ) 



(24) 



it is possible to obtain an asymptotic estimate also for 
E corr (r s — > oo). One can perform the same compar- 
ison using the LRDMC data. In Fig[3] the value of 
Ecorr/ E ex is reported as function of r s for all b stud- 
ied (b = 0.1,0.3,0.5,0.75,1,2,4). The points have been 
fitted using the following function: 



Ecorr s) 

E ex (r s ) 



ln(r s ) 



c 



(25) 



where a, b, and c are fitting parameters. According to 
Eq. [Ml the slowest decaying contribution to E C orr/E ex 
should be oc l/ln(r s ), which has been included in the 
fit. Therefore an accurate extrapolation to r s — > oo is 
very difficult, since the logarithmic correction in the low 
density regime is extremely slow, and in order to have a 
reliable estimate of a, one should compute the correlation 
energy at much higher values of r s . For instance, in the 
3SRA-STLS theory^, the ratio E corr /E ex converges to 
the value in Eq. [221 for r s > 1000, which is of course 
a regime where the QMC framework is not ergodic. In 
any case, in the range of densities from r s = 15 to r s = 
50 (the lowest density taken into account) the function 
in Eq. [5S] fits well our data points, apart from the case 
with b = 4, which requires even lower densities to enter 
into its asymptotic regime. In our fits, 5, assumes values 
ranging from 1.0 to 1.1, which are slightly larger than 
those yielded by the STLS theory. Thus also our data 
seem to support the idea that the correlation energy in 
quasi-one dimensional systems behaves as 



E CO rr (r s — > OO) OC - 



hi(r s ) 



(26) 



or at least this behavior is compatible with our data until 
r s = 50. 

A good parametric representation of the correlation 
energy has to satisfy the asymptotic behavior both at 
r s = and r s = oo. Agosti et used the interpolation 
formula 



e(r s 



C + Dri 



(27) 



with p, q, C , and D variational parameters, in order to fit 
their STLS data for a similar 1DEG model system (with 
hard confinement of a 2DEG in one direction). They 
found the above function fits correlation energies very 
well in the intermediate range of densities, but it is ap- 
parent that their expression does not yield the "exact" 
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FIG. 3: (Color online) Values of E CO rr/E ex plotted versus 
r s in the low density regime for various thickness b. The 
curves are obtained from a fit of the LRDMC points using 
the function in Eq. [25] in the range 15 < r s < 50. 



asymptotic behavior at high and low densities. Here we 
give a simple analytic representation of the correlation 
energy which fits our simulation data over a wide range 
of densities (from r s = 0.05 to r s = 50) and for different 
values of lateral confinement. Moreover it includes the 
correct low and high density limits. Following the same 
lines as Perdew-Wang^ 4 - and Attaccalite et ai4^, we use 
a parametric function which reads 



A + Br 7 } + Ck 



\n(l + ar s + fir?). (28) 



The total number of parameters is 7: 3 linear coefficients 
(A, B, and C) and one exponent (0 < n < 2) for the 
polynomial part, 2 linear coefficients (a, and (3) and one 
exponent (m > 1) in the argument of the logarithm. The 
high density limit of the parametrization is 



e(r s -» 0) 
and its low density limit is 



a 
A 



m ln(r s ) 
e(r s ^oc) = -- — — . 



(29) 



(30) 



Thus the expected behavior of the correlation energy is 
correctly reproduced by our parametric form (Eq. I28[) for 
both the high and low density limits. In order to obtain a 
priori the exact high density result (Eq. I22p known from 
the RPA theory (see Appendix [A"]) , we fix the ratio a/ A 
to be equal to A corr /(ir 4 b 2 ). Therefore the number of 
independent parameters in the function (|28p is reduced 
to 6. In particular, we determine A from the high density 
limit, while the other parameters are free to minimize the 
X 2 - Their optimal values are listed in Tab. HVl 

We define the accuracy r\ of the parametrization by 



1 M 

8=1 



.(*)-£(*)!, 



(31) 



where M is the total number of r s points computed for a 
given b. In practice, r\ is the average of the residuals and 
measures the discrepancy between the computed value 
E corr and its parametric value e. The order of magnitude 
of the accuracy r\ is between 10 -4 Ryd* and 10 -5 Ryd* , 
depending on the thickness of the wire. It means that the 
parametrization is very accurate in the range of density 
with < r s < 50 and for values 0.1 < b < 4 of the 
width parameter. In Fig. [4] we plot the points and the 
parametrization curves for the correlation energy in the 
range of density and width taken into account. Notice 
that the thiner is the wire, the more correlated is its 
ground state. 
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FIG. 4: (Color online) Correlation energy versus the den- 
sity parameter r s . The points are the LRDMC data used 
in the fit. Their dimension is much bigger than the er- 
ror bars. We computed the LRDMC correlation energy 
at r s = 0.1, 0.2, 0.3, 0.4, 0.6, 0.8, 1, 2, 3, 4, 6, 8, 10, 15, 20, 35, 50. 
The curves are the parametrization of the correlation energy 
for various values of the width parameter b. In the inset, the 
semi-log plot magnifies the region around the minimum of the 
energy correlation functional. 

To conclude this section, in Fig. [5] we compare our 
LRDMC results with those obtained by Calmels and 
Gold using the 3SRA approach^. The relative difference 
between the two methods is about 20% around the mini- 
mum, where the absolute value of the correlation energy 
is larger. The STLS method with the 3SRA sum rule is 
able to yield a large fraction (80% or more) of the cor- 
relation energy. However our results, which are formally 
exact, represent a further improvement with respect to 
the previous estimate of the correlation energy for this 
model quantum wire. 



VI. PAIR CORRELATIONS 

In this section we study the density-density and spin- 
spin correlation functions, which are useful quantities for 
assessing the nature of the ground state of the system. 
We first analyze the charge g PP (r) and spin g^^ir) pair 
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TABLE IV: Optimal fit parameters for the correlation energy, as parametrized in Eq. 1281 Different values of b are taken into 
account, for each of them we give its parametrization. In the last rows, we report the reduced x 2 and the overall numerical 
accuracy r\ in Ryd* , defined in Eq.|3l] A has been determined from the high density limit 1221 which is therefore exactly fulfilled 
by our fit. The error of the parameters in the last digits is reported in parenthesis. 





b = 0.1 


b = 0.3 


b = 0.5 


b = 0.75 


b = 1.0 


b = 2.0 


b = 4.0 


A 


4.66(5) 


9.50(9) 


16.40(11) 


22.53(24) 


32.1(1.3) 


110.5(2.8) 


413.0(5.9) 


B 


2.092(24) 


1.85(17) 


2.90(13) 


2.09(12) 


3.77(23) 


7.90(49) 


10.8(7) 


C 


3.735(34) 


5.64(8) 


6.235(46) 


7.363(44) 


7.576(36) 


8.37(6) 


7.99(20) 


n 


1.379(10) 


0.882(13) 


0.908(7) 


0.906(8) 


0.941(49) 


1.287(26) 


1.549(35) 


a. 


23.63(26) 


5.346(53) 


3.323(22) 


2.029(22) 


1.63(7) 


1.399(35) 


1.308(19) 




109.9(5.5) 


6.69(43) 


2.23(7) 


0.394(12) 


0.198(8) 


0.0481(19) 


0.0120(9) 


m 


1.837(18) 


3.110(42) 


3.368(23) 


4.070(30) 


4.086(21) 


4.260(25) 


4.165(40) 
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FIG. 5: (Color online) In the lower panel, we plot the corre- 
lation energy obtained by Calmels and Gold^ (points) and 
the parametrization of our LRDMC data (curves) for various 
values of the width parameter b. In the upper panel, we com- 
pute the relative difference |e — Ec a imeis\l 'e interpolated all 
over the density range taken into account. 



correlations, defined as 

g pp (r) = 1/2 (g n (r)+g n (r)), (32) 
9aa{r) = l/2(g n (r)-g n (r)), 

where we used the spin resolved pair distribution function 

rf-0), (33) 



with a, (3 =T,|, and — Pi — p/2 is the density of 
the two spin components in the unpolarized system. We 
computed g PP {r) and gaaij^ of the wires with thickness 
b = 0.1 and densities r s — 1,2,4, by means of the 
LRDMC algorithm and the forward walking technique^, 
which yields unbiased expectation values on the ground 
state of the system also for operators which do not com- 
mute with the Hamiltonian. The correlation functions 
are drawn in Figs. and [7] 

The charge-charge correlation functions in Fig. [S] re- 
veal the strong effect of the electronic correlation in the 
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FIG. 6: (Color online) Charge-charge pair correlation function 
g PP (r) computed from LRDMC simulations with b = 0.1 and 
N = 82 at different densities. In the inset, the short-range 
part of g PP (r) is magnified. 
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FIG. 7: (Color online) Spin-spin pair correlation function 
gao(r) computed from LRDMC simulations with b = 0.1 and 
N = 82. 



low density regime. As the density is decreased, the par- 
ticles repel each other with an effective interaction which 
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is stronger. Consequently g PP (0) is smaller, as one can 
see in the inset of the Fig. [Sj At the same time the fluc- 
tuations in the g PP {r) are larger, with periodicity 2r s OQ, 
and a slow decay. We will see however that this decay 
is not slow enough to give rise to a true long range or- 
der (Wigner crystal) , since the quantum fluctuations will 
prevent the freezing of the charge in ID. On the other 
hand, the short-range behavior of the spin-spin correla- 
tion functions in Fig.[7]shows the antiferromagnetic char- 
acter of the coupling among electrons, and the underlying 
periodicity in the spin sector is 4r s ag, i.e. twice the mean 
interparticle distance. 

The analytical behavior of charge and spin correlations 
for the quasi 1 DEG with long range interaction has been 
obtained by SchulzAS using bosonization techniques ap- 
plied to an effective one dimensional Hamiltonian with 
linearized kinetic energy. For that model Hamiltonian, it 
turns out that the charge correlation function exhibits a 
slow decay of its Akp component 

(p(x)p(O)) ~ Acos(4:k F x) exp(-4<Vhi(x)), (34) 

with c an interaction dependent parameter. Its behavior 
has been related to a quasi order of the electrons, since 
in one dimension there is no true long-range order. Their 
fluctuations are almost frozen to maximize the interpar- 
ticle distance along the wire, and their relative positions 
are pinned around lattice sites with periodicity 2r s ag, ex- 
actly as we have seen in the g<ja(r) for our model. The 
formation of this quasi Wigner crystal comes from the 
1 jr tail of the potential, since in the case of short-range 
interactions the Akp component of the charge correlation 
function decays much faster. On the other hand the spin 
correlation function has no singularity at Akp and ex- 
hibits the slowest decay for the 2k p component. Indeed, 
according to the bosonization technique^, the large dis- 
tance spin correlations are given by: 

(a(x)a(Q)) ~ B cos(2kpx) exp(— cy/ln(x))/x, (35) 

where c is the same as in Eq. 1341 

In order to check whether the Eqs. [34l and [35l are valid 
for our quantum wire model with long-range interactions 
and quadratic dispersion, we have also studied the charge 
S pp (k) and spin S aa -(k) structure factors, defined as 

S pp (k) = {p{k)p(-k))/N (36) 
S aa {k) = (a(k)a(-k))/N, 

where p(k) (cr(k)) is the Fourier component of the local 
charge (spin) density. We have computed the structure 
factors at several k values for b = 0.1, r s = 1,2, 4, using 
N = 10,22,42,62,82 in the LRDMC calculations and 
N = 10, 22, 42, 82, 162 in the less expensive VMC calcu- 
lations. In Figs. [5] and [5] we plot S pp (k) and S aa (k) for 
the LRDMC simulations with the largest number of par- 
ticles (N = 82). While the spin-spin correlations show a 
peak at k = 2k p, for the charge degrees of freedom the 
highest peak arises at k = Akp, which corresponds to the 
periodicity of the quasi Wigner crystal. 
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FIG. 8: (Color online) Charge structure factor S pp (k) com- 
puted from LRDMC simulations with b — 0.1 and N = 82. 
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FIG. 9: (Color online) Spin structure factor Sa-a{k) computed 
from LRDMC simulations with b = 0.1 and N = 82. 

In particular, we have studied the dependence of the 
peak heights S pp (4:kp, N) and S aa (2kp, N) on the num- 
ber of particles N, and compared it with the behavior 
predicted by the bosonization. From Eqs. [34] and [35l 
one can easily obtain^: 

S pp (Ak F ,N) = a x L exp(~4c ^log L) + a 2 , (37) 
S aa {2k F ,N) = a 3 (yioiL/c+l/c 2 )exp(-cV!oiI) 
+ a 4 , (38) 

with L = 2r s N , 01,02,03,04 model and density dependent 
parameters, and c the same as in the bosonization results. 
We have then used Eqs. [38l and [37l to fit our results, ob- 
taining first ai, 02, c from the fit of S pp (Akp, N) and then 
03, 04 from Scro-(2kp, N), with c fixed by S pp (AkF, N). In 
Figs. [TU] and QT] we plot the curves which best interpo- 
late the VMC and LRDMC values for S pp {Ak F ,N) and 
S aa (2k F ,N). 

Our VMC and LRDMC results seem consistent with 
the predictions of the bosonization, at least in the range 
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FIG. 10: Peak of the charge structure factor at 4fcp versus the 
system size for r s = 1, 2, 4, b = 0.1, obtained from VMC and 
LRDMC calculations. The curves fit the points with function 
in Eq.rjT] 



CO 




20 40 60 80 100 120 140 160 
Number of particles 



FIG. 11: Peak of the spin structure factor at 2Uf versus the 
system size for r s = 1, 2, 4, 6 = 0.1, obtained from VMC and 
LRDMC calculations. The curves fit the points with function 
in Eg. [381 



of system sizes taken into account. As it is apparent from 
Figs. [10] and [TTJ the points are well fitted by Eqs. [37] 
and [35] respectively, the distance from the interpolating 
curves being usually less than 2 standard deviations. As 
already reported in Ref. H3, VMC results are in quali- 
tative agreement with the bosonization findings, and the 
variational ansatz in Eq.[l0]is good enough to capture the 
correct ground state properties. Indeed the LRDMC pro- 
jection changes only quantitatively the VMC points, by 
reducing the charge structure factor, and by enhancing 
the spin structure factor, which however remains finite in 
the thermodynamic limit in accordance with Eq. 1381 We 
emphasize once again that for this one dimensional sys- 
tem the LRDMC results are "exact" within their statis- 
tical accuracy, since the LRDMC yields the exact ground 
state energy and the points in FigfTOlandfTTlare obtained 



using the forward walking technique^, which provides an 
unbiased expectation value for each correlation function. 

The behavior of the spin and charge structure factor at 
small momenta reveals important features of the ground 
state, which are related to the low energy modes of the 
system. According to the bosonization results^, S pp (k) 
should behave as oc \k\/ */] ln(fc)|, while S crcr (k) should go 
linearly with Our results are plotted in Figs. [T2"land 
1131 In particular, in Fig. [12] we draw the renormalized 
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FIG. 12: (Color online) Plot of f p (k) 
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versus fc/fci? for r s = 1,2,4, b = 0.1, obtained from VMC 
(N = 162) and DMC (N = 82) calculations. The horizontal 
line is drawn as an eye guide. 



charge structure factor f p (k) = -4= 



2 V|ln(fc)| 



S pp (k). Both 



the VMC and LRDMC data are in agreement with the 
small k limit of S pp {k): 



lim S p p(k) = a g 

k— *0 



|fej 



V\HW\' 



(39) 



where a g is a factor which seems close to 1 and only 
slightly dependent on r s , although the logarithmic be- 
havior of f p (k) would require much smaller k in Fig. 1 121 in 
order to have an accurate extrapolation for a g . Anyway, 
our result agrees at least qualitatively with the bosoniza- 
tion findings. Moreover, it is possible to derive ana- 
lytically the behavior in Eq. [35] directly from our vari- 
ational wave function ^t, using the expression obtained 
by Reatto and Chester—, that relates the small momenta 
behavior of S pp (k) with the 2-body Jastrow factor in- 
cluded in ^t- 



S pp (k) 



S Q (k) 



1 + 2pu(k)S {k) ' 



(40) 



where u(k) is the Fourier transform of the Jastrow func- 
tion. The above relation is approximate for finite fc, but 
it becomes exact to the leading order in fc, for k — > 0, for 
the variational structure factor, i.e., for the one evaluated 
from the Slater-Jastrow trial function. After plugging the 
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definition of u(k) (Eqs. [Til and [T2|) into the above equa- 
tion, in the limit k —> we have 



S 



m 1 ^ i fc l 
"" W 7 2 ^TnW' 



(41) 



where 7 is the parameter of the optimized RPA Jastrow 
factor. It is therefore clear that our choice of u(k) sat- 
isfies the correct behavior of the charge structure factor 
already at the variational level, and the LRDMC simula- 
tions do not change this behavior (see Fig. [T2"]), Moreover, 
from the knowledge of S pp (k) we can infer the behavior 
of the low energy charge excitations (plasmons) oj p {k). 
Indeed, a variational estimate of u> p (k) is given by (see 
Appendix [B| 



Wp(fc) 



S pp (k) 



and in the limit A; — ^ it turns out that 

2 



>(*0 = 



\k\y/\h^k\ 



(42) 



(43) 



This expression should be compared with the low energy 
spectrum provided by bosonization studie o 16 i ' 50 of the 
Coulomb Luttinger liquid, which reads 

ui p {k) = v p (k)\k\, 



Vp(k) = W(l + - 9i + 2V(k)/irv F ), (44) 

where v p is the charge velocity, and g\ is the amplitude of 
the backward scattering process. For small k excitations, 
we have 



and it is the same behavior as the spin structure fac- 
tor of a non-interacting gas. Therefore, the interaction 
leaves unchanged the S aa (k) tail at small k, which is 
another striking feature of the spin-charge separation. 
Notice that if we use S aer (k) to estimate the low lying 
spin- wave excitations (see Appendix iBj). we will obtain a 
spin velocity equal to Vp, and independent of r s . This is 
of course a variational estimate, since the true spin veloc- 
ity strongly depends on the density and is significantly 
reduced by the effective interaction^, being equal to vf 
only in the high density weak interaction regime. 
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FIG. 13: (Color online) Plot of f a {k) = S tTCT (fe) versus 
k/k F for r s = 1,2,4, 6 = 0.1, obtained from VMC (N = 162) 
and DMC (N = 82) calculations. The horizontal line is drawn 
as an eye guide. 



limw p (fc) = VT+^T-— \k\^/\h^k\, (45) 

k— »0 



which corresponds to our findings in Eq. 1431 Since in one 
dimension the long-wavelength spin and charge modes 
are independent, they have different velocities. This dif- 
ference, due to the so-called spin-charge separation, has 
been seen in a remarkable experiment by Auslaender et 
alA, and predicted by the Luttinger liquid theory. In- 
deed, according to this theory the spin excitations are 



LUa(k) = vfJI -9i\k\- 



(46) 



The spin dispersion is linear, since it is not affected by 
the long-range tail of the Coulomb interaction, in con- 
trast to the charge velocity v p , which is renormalized by 
the Fourier transform of the potential at small k. The 
linear behavior of the spin branch is reflected by the lin- 
ear decay of S aa {k) as k goes to 0. In Fig. [13] we plot 
the renormalized spin structure factor f a (k) = !1 k E -S aa (k) 
computed by carrying out both VMC and LRDMC sim- 
ulations for r s — 1,2,4. Since the value of f a (k) is 1 in 
the limit of small fc, it turns out that 



Scrcr (k) 



k 

vf ' 



(47) 



VII. CONCLUSIONS 

In this paper we have carried out extensive Monte 
Carlo simulations to compute the ground state properties 
of the quantum wire model with unscreened long range 
interactions. We have used the novel LRDMC frame- 
work, which is shown to be more efficient than the stan- 
dard DMC algorithm in the strong coupling regime, i.e. 
at low densities, when the exchange is extremely small 
and the features of a quasi- Wigner crystal are manifest. 
We computed the exact correlation energy, and found a 
simple and accurate parametrization, which fits the cor- 
relation energy over a wide range of electron densities 
and lateral confinements. This parametrization includes 
the correct behavior at high densities (e(r s — > 0) oc r 2 s ), 
given by the RPA approximation. On the other hand, 
we guessed the asymptotic behavior at low densities from 
our data, and we found that e(r s — > 00) cx — ln(r s )/r s fits 
well the LRDMC correlation energies in the strong cou- 
pling regime. We believe that our parametrization pro- 
vides an extremely reliable functional for further DFT 
computations of quasi one dimensional systems. Last 
but not least, we showed that the pair correlations of 
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our model, exactly computed by means of LRDMC sim- 
ulations and forward walking techniques, reproduce all 
features of the so-called Coulomb Luttinger liquid, i.e. 
a Luttinger liquid with linear dispersion and long-range 
Coulomb interaction. In particular, our results are com- 
patible with the slow decay (oc exp(— 4c-y/hi(x)) of the 
spin structure factor at 4fcp, which is the signature of a 
quasi order of the charge degrees of freedom. Moreover, 
the small fc behavior of the static structure factor is in 
agreement with the bosonization findings, both for the 
charge and spin modes. 

We plan to extend this study also to the spin polarized 
case, and to provide a spin resolved energy functional to 
be used as input of one dimensional DFT calculations 
with local spin density approximation. 
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APPENDIX A: RPA CALCULATION OF THE 
CORRELATION ENERGY 



and the physical density-density response function of the 
free ID electron gas is 

lim x°(k,uj + iri). (A3) 

Following the seminal work of Gcll-Mann and 
Brueckner— , whose approach has also been used 
by Rajagopal and Kimball 53 for the 2D case, we define 
the so called electron-hole propagator 

dk / dt f(k)(l- f(k + q)) 

-co J — CO 

e-^eM~\t\(\q 2 + kq)), (A4) 
where f(x) — 9(\x\ — 1) is the zero temperature Fermi 
distribution, with 9 the step function: 

It is easy to see that the response function in Eq. lA2l can 
be simply related to the electron-hole propagator by the 
relation 



In this appendix, we compute the correlation energy in 
the high density limit using the random phase approxi- 
mation (RPA) . The RPA correlation energy reads 5 -^ 



E. 



RPA 



L 

2^ 



+ 00 



dk E(k), 



1 Ifcl f +0 ° 

E ( k ) = T-y / dx Mi-v(kb) x °(k,ik\)) 



4tt N J 

+v(kb)x°(k,ik\), (Al) 
where v(kb) is the Fourier transform of the potential, 



1 



X u (fc, W ) = — In 



lo 2 - (fc 2 - v F kf 



2nk \lu 2 - (fc 2 + VFky 



(A2) 



2-Kk, 



-Qq(u), 



(A6) 



where q and u are dimensionless variables defined by 



fc = fci?g 

lo = i kpq vfu. 



(A7) 



Therefore, if we rewrite Eg. lAll in terms of the electron- 
hole propagator and using the dimensionless variables u 
and q, we obtain 



+ °° r+ °° r ' ar s ( qb\ „ , \ ar s ( qb 

au q < m i i -t- ' ' 

— \" n( 'I 1 ' 



= 2^ Jo ^ Jo ^ q I" 1 { 1 + ^ V {^) Q * {U) )-^ V & QM 



Y /"TOO P~\-00 -^'^j 

—— rij / dq I du q) ( , 

27r{ar s ) z J Q J ^ n \ 2n J \ar. 



?(«), (A8) 



where a = A/n in ID. To pass from the first line to 
the second one we have Taylor expanded the logarithm. 
Notice that in ID, in contrast to the 3D and 2D cases, 



the integrals converge in all orders of the expansion, since 
v(qb) (= Vb(q) in Eq. [3]) diverges only logarithmically at 
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small q: 



v{qb) 



-41n(g6) if q 



Moreover, the leading order in the r s expansion for 
Eq. IA8I is given by n = 2 (i.e. the direct lowest ring di- 
agram in the perturbative expansion of the interaction). 
Since we are interested in the lowest order r s expansion 
of the RPA correlation energy, we keep the term with 
n = 2 and we discard the others. Thus, we get 



where we used the asymptotic behavior of the poten- 
tial v(x) for x — > +oo (Eq. IA9j) . and defined A exc h = 

(A9) Jo °° dz z v(z). It is apparent that the exchange sec- 
ond order diagram contributes only to the fourth order 
of r«. Therefore, at the lowest order in r s , the correlation 



energy is: 



A 
'tt 4 6 2 ' 



(A16) 



E. 



RPA 



\j^r dqqv2 {t) F{qi (Mo) 



with F(q) = f °°du Qq(u) the integral over the dimen- 
sionless frequency of the 2-particles electron-hole propa- 
gator at a given momentum transfer q. Indeed F(q) can 
be explicitly written as 



f(q) = — [ dh [ dk 2 

1 Jl-q Jl-q 



1 



q(ki +k 2 )' 



(All) 



and for zero q-transfer F(Q) = ir. If we rescale the vari- 
able q in the integration IA10I (q — > ^jf-q), and we keep 
the lowest order in r s , the RPA correlation energy reads: 



with A = Jq dz z v 2 (z) = 4.9348. This result turns 
out to be the same as a high density extrapolation^ 2 - of 
the correlation energy for the same model studied here, 
obtained by Gold and Calmels within the so called mean 
spherical approximation (MSA) . In two and three dimen- 
sions, the MSA yielded high density expressions of the 
correlation energy which were slightly different from the 
RPA findings 23 . In this case however, the MSA and the 
RPA results are in perfect agreement. 



APPENDIX B: VARIATIONAL ENERGIES OF 
CHARGE EXCITATIONS 

Let \^o) be the ground state of the Hamiltonian 



E. 



RPA 



(47T 



ar & 
,2 [ — 



+ 00 



dzzv\z). (A12) ff = ^e(fe)4 i<TCfc , ff + i^y( g )(p_ 9 p ? -iV), (Bl) 



k.a 



In order to make sure that this is the correct high energy 
limit of the correlation energy, we need to consider also 
the second order exchange contribution in the perturba- 
tion theory. It is neglected in the RPA approximation, 
but can yield non trivial corrections in the r s expansion, 
like in the two^ 3 - and thread dimensional electron gas. 
The second order exchange is 

KL h = \ l +C ° dq qv (£) F exc ,(q), (A13) 
where now 



where e(k) is the dispersion of the non interacting sys- 
tem, V(q) is the Fourier transform of the interaction, 



and p q — J2k a c k+q a Ck ^ ^ s ^ ne Fourier transform of the 
charge density operator. In analogy with the Feynman's 
construction for the liquid Helium 54 , a variational wave 
function for a charge excitation (plasmon) with momen- 
tum q is given by 



P?l*o) 



(B2) 



Fexch{q) = — / dki I dk 2 

1 Jl-q Jl-q 



Its variational energy (& q \H |* g )/(^ ? |* g ) is E q , while 
the GS energy is Eq. Notice that the normalization of 



l* 9 > is 



q 2 + q(ki + k 2 



A14) 



(q + ki +k 2 )b 



One can easily see that F exc h{0) = tt v The dif- 

ference from the second order direct ring (Eqs. IA10I and 
lAlip is the vertex interaction V{q) computed at q+ki+k 2 
instead of q, and the overall factor is reduced by a factor 
of 2, due to the spin summation which is now restricted 
only to the parallel contribution. If we rescale the vari- 
able q as before (q — > ^jf-q), we find 



(*,!*,) = NS pp (q)(* \*o), 



(B3) 



where S pp (q) = jj(p- q p q ) is the static charge structure 
factor. We are now going to find an expression which 
allows us to estimate the excitation energy of a plasmon 
with momentum q from knowledge of S pp (q). We start 
by evaluating the double commutator 



[[Pq,H},p-q] = ^{e k+ q+e k -q - ^ k )c{ 



k Ck, 



(B4) 



E. 



ii 



1 1 



2 (4tt) 2 



dz z v(z) v ( 



where we used the Hamiltonian in Eq. IBll On the other 
hand, we have 



4A 



exch 4 
T 



7T 6 6 4 



(A15) 



(Wo|[[p g ,g],p- g ]|^ ) 

(*o|*o) 



2NS pp {q){E q -Eo), (B5) 
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by applying the definition of the double commutator 
[[Pq, H],p- q ] = p q Hp_ q + p- q Hp q - p q p- q H - Hp q p_ q . 
Merging Eq. IB4I and Eq. IB5[ we are led to the following 
identity: 



E„ — En 



2NS pp (q){* \* ) 



(B6) 



Eq — Eq is an estimate of the plasmon excitation with 
momentum q. In general the ansatz \^ q ) = Pql^o) is not 
exact for the lowest energy wave function with momen- 
tum q, but gives a variational energy, since it belongs to 
the same q-subspace as the true excited state and is or- 
thogonal to subspaces with different q' . In the limit of q 
small, Eq. IB6I turns out to be 



E q - E 



(*0|£*3fo4<*l*0) q 2 



2JV(* |*o) 



S pp (q) 



where we used the quadratic dispersion e(fc) = k 2 in Ryd* 
units. Therefore, from the knowledge of S pp (k) calcu- 
lated for the GS of the system we can evaluate its excita- 
tion spectrum in a variational way. Moreover the smaller 
q is, the better Eq. IB7I approximates the true plasmon 
energy. In the same way, one can estimate the energy 
of the spin- wave excitations (spinons), by using the re- 
lation in Eq. IB7I with S pp {k) replaced by the static spin 
structure factor S aa {k). 



Spp(q)' 



(B7) 
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